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We study a generic class of inelastic soft sphere models with a binary collision rate g" that 
depends on the relative velocity g. This includes previously studied inelastic hard spheres (y = 1) 
and inelastic Maxwell molecules (y = 0). We develop a new asymptotic method for analyzing large 
deviations from Gaussian behavior for the velocity distribution function /(c). The framework is that 
of the spatially uniform nonlinear Boltzmann equation and special emphasis is put on the situation 
where the system is driven by white noise. Depending on the value of exponent v, three different 
situations are reported. For v < —2, the non-equilibrium steady state is a repelling fixed point of 
the dynamics. For v > —2, it becomes an attractive fixed point, with velocity distributions /(c) 
having stretched exponential behavior at large c. The corresponding dominant behavior of /(c) is 
computed together with sub-leading corrections. In the marginally stable case v — —2, the high 
energy tail of /(c) is of power law type and the associated exponents are calculated. Our analytical 
predictions are confronted with Monte Carlo simulations, with a remarkably good agreement. 
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I. INTRODUCTION 

The interest in kinetic theory of dissipative systems, such as granular gases and fluids Q, 0, 0, 0, has caused a 
great revival in the study of the Boltzmann equation 0, Q . Not surprisingly, the introduction of Maxwell models with 
their energy independent collision rate -which simplifies the nonlinear collision term to a convolution product- has 
had a great part in that H,0>E3 in this paper the focus is on the velocity distribution function (v.d.f.) F(v,t) in 
spatially uniform states of inelastic systems, evolving according to inelastic generalizations of the Boltzmann equation 
for classical repulsive power law interactions. For these systems we study the asymptotic properties of the v.d.f.'s 
at large times and at large velocities. This will be done for cases without energy supply, i.e. freely cooling systems, 
as well as for driven systems. The latter ones may approach a non-equilibrium steady state (NESS), and the former 
ones approach scaling states, described by scaling or similarity solutions of the nonlinear Boltzmann equation. Both 
types of asymptotic states show features of universality, such as independence of initial states, and independence of 
the strength of the energy input, but do depend on the type of driving device. 

The big boom occurred in 2002 after the discovery of an exact scaling solution, /(c) = (2/7r)(l + c 2 )~ 2 of the 
spatially uniform nonlinear Boltzmann equation (NLBE) for a freely cooling one-dimensional Maxwell model Il2| . 
Here, the velocity distribution function (v.d.f) has the scaling form F(v,t) = Vn d f(v/yn) where vo(t) is the r.m.s 
velocity and d is the number of spatial dimensions. Subsequent analysis |13l ll4l Tl5l fl q | has shown that the NLBE 
for freely cooling inelastic Maxwell models in d dimensions has a scaling solution with a power law tail /(c) ~ l/c s . 
The power law exponent s can be calculated from a transcendental equation, and depends on the dimensionality d 
of the system and on the coefficient of restitution a, where 1 — a 2 measures the fractional energy loss in an inelastic 
collision. The proof that the v.d.f. F(v, t) for general initial values F(v, 0) approaches for large t this scaling solution 
has been given in ^tJ- The scaling and NESS solutions mentioned above have the remarkable property of possessing 
a power law tail, /(c) ~ l/c s , which is highly overpopulated at large velocities as compared to a Maxwellian v.d.f 
~ exp(— c 2 ). Solutions with overpopulated tails of stretched exponential form, /(c) ~ exp(— /3c h ) with < b < 2, have 
been studied before both analytically and by numerical simulations (both Molecular Dynamics and Monte Carlo) for 
inelastic hard spheres jEEEIllGSIMIMIIllIllMIll, as well as for inelastic soft spheres 0, |2(| . Such states 
appear not only in freely cooling systems but also in driven systems. These systems are in principle more interesting 
because the collisional loss of energy may be compensated by spatially uniform heating devices. They may drive the 
system into a NESS, with possibilities of experimental verification. 

Many different types of heating devices have been invented: 

i) white noise, which adds random velocity increments to the particles HHmQHEHIHIillii] 
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ii) a deterministic friction force a = 7ov|u| e with 9 > and v = v/v (nonlinear friction model). For 9 = 1 it gives 
the Gaussian thermostat equivalent to a linear rescaling of the velocities, and for 6 = it gives the 'gravity' 
thermostat, which models something like gliding friction These thermostats have been used in theoretical, 
molecular dynamics and Monte Carlo studies fla Hij. Ell l22l I2H I30L l3l | to add (70 > 0) or to remove (70 < 0) 
energy from the system. 

Hi) a heating device, recently proposed in Ref. |32| . that feeds energy into the system in the ultra high energy tail 
of the v.d.f and forces the system into a steady state with a power law tail. We shall come back to this model 
at the end of section IVll 

iv) a heat bath consisting of an ideal gas kept in a thermal equilibrium state |33| . or a more exotic device that 
maintains the bath in some non-equilibrium steady state in which the v.d.f. is, for instance, described by a Levy 
distribution, having power law tails |34|. 

Over-populated tails, corresponding to large deviations from the Maxwellian v.d.f, have also been observed in many 
experimental studies of driven granular gases and fluids |35j . The heating devices used there are difficult to model in 
a simple kinetic theory context 36] . 

The existence of such tails is remarkable because the standard NLBE for classical gases and fluids with conservative 
interactions [37LI38L l39t l40| provides the basic notion of rapid relaxation of a general initial distribution F(v, 0) towards 
a steady state described by a Maxwellian v.d.f. Here the steady state solutions satisfy the detailed balance relation 
F(v[, oo)F(v' 2 , 00) = F(vi, oo)F(v 2 , 00) between pre-collision velocities (v l5 v 2 ) and post-collision ones (v'^vf,). The 
detailed balance relation in combination with the if -theorem forces the steady state solution to be a Maxwellian. 
The basic reason for the large deviations from Maxwellian behavior is the violation of detailed balance in dissipative 
collisions, together with the breakdown of the ii-theorem |15|. 

The goal of this paper is to develop a new asymptotic method for analyzing the large deviations from Maxwellian 
behavior in the high energy tails of the v.d.f. observed in NESS solutions of the nonlinear Boltzmann equation for 
dissipative systems, driven by an energy source. The class of v models studied are the inelastic soft spheres prj I32] ] 
with a binary collision rate, scaling like g u , where g is the relative speed. It includes the previously studied inelastic 
hard spheres (v = 1) and inelastic Maxwell molecules (y = 0). 

The method can be applied to all driving devices listed above. Depending on the device used, the models have in 
general either i) a stable NESS, i.e attracting fixed point solution, for b > or v > v c with stretched exponential tails 
/(c) ~ exp(— (3c b ) and higher asymptotic corrections, like /(c) ~ c x exp(— (3c b + (3'c b ), ii) a marginally stable NESS 
for a threshold model with b — or v = v c with power law tails, and Hi) an unstable NESS, i.e. a repelling fixed point 
solution for b < or v < v c . For free cooling and Gaussian thermostats b — v, for white noise driving b = 1 + \v and 
for nonlinear friction models b = v + 1 — 6. The preliminary results have been published in Ref. 41] . The present 
paper focusses on the mathematical method, and on the simple application of white noise driving (i) , and on the ultra 
high energy source (iii). The remaining results will be published elsewhere | l2j. 

The crux of the new asymptotic analysis is the construction in Section [H] of a linearized collision operator, whose 
eigenfunctions are powers of the velocity c, and whose eigenvalues determine the power law exponent in the tail of 
the distributions. The Fourier transform method used in Refs. [l4l Il5j for determining power law tails can only be 
applied to Maxwell models (y — 0). In Section ITTT1 we study the spectral properties of the aforementioned operator, 
with details given in appendices El and These two sections plus appendices can be considered as the generalization 
to inelastic soft spheres of the mathematical theory for inelastic Maxwell molecules in Refs. 0,@. In Section llVl we 
derive the energy balance equation for the white noise driven case, determine the stability of the fixed point solutions, 
and derive the integral equation for the scaling form /(c). In Section^we present the high energy tails of exponential 
type, and in Section IVII those of power law type, where also the ultra high energy source is discussed. The results 
are supported by Direct Monte Carlo Simulations (DSMC scheme 0]). We end in Section IVTT1 with conclusions and 
perspectives. 

II. BASICS OF INELASTIC SCATTERING MODELS 

The nonlinear Boltzmann equation for dissipative interactions in a spatially freely evolving state can be put in 
a broader perspective, that covers both elastic and inelastic collisions, as well as interactions where the scattering 
of particles is described by deterministic (conservative or dissipative) forces or by stochastic ones. To do so it is 
convenient to interpret the Boltzmann equation as a stochastic process, similar to the presentations in the classical 
articles of Waldmann |44|. Uhlenbeck and Ford 0], or m Ulam's stochastic model |4^. The last one shows the basics 
of the approach of a one-dimensional gas of elastic particles towards a Maxwellian distribution. 

We consider a spatially homogeneous fluid or gas of elastic or inelastic particles in d dimensions, described by 
an isotropic velocity distribution function, F(y,t) — _F(|v|,£), and evolving according to the nonlinear Boltzmann 
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equation, 

d t F(v, t) = I(v\F) = J dwdv'dw' J dn [W(v, w|v', w'; n) 

xF(v', t)F(w',t) - W(v', w'|v, w; n)F(v, t)F(w, t)} , (2.1) 

where the binary collisions are described through a transition probability per unit time, W(v', w'|v, w; n). The 
loss term 7i OS s and the gain term / ga in represent respectively the contributions from the direct collisions (v, w) — > 
(v',w'), and from the restituting collisions, (v',w') — > (v, w). Here the direct and restituting velocities have been 
parameterized in terms of the incoming velocities (v, w), and an impact (unit) vector n, that is chosen on the surface 
of a unit sphere with a probability proportional to the collision frequency K(g±, gu). It depends only on the length 
g± = |gj_| and |<7||| of the vector arguments. We further use the notations 0|| = a ■ n and = a— a||n, and a = a/a 
is a unit vector. 

The transition probability for the inelastic collisions, (v, w) — > (v', w'), is given by (2(|, 

M/(v',w'|v,w;n) = K(g ± , g^^G' - G)S^\^ ± - + <*g\\), (2-2) 

where G = i(v + w), g = v — w, and a obeys < a < 1. In the subsequent analysis, we refrain from explicitly 
indicating the dimension of the argument for the Dirac functions. The value a = 1 corresponds to the elastic case, 
and a = to the totally inelastic case. The elastic case, possibly driven by by source(s) and/or sink(s), will not be 
considered here. On the other hand, the quasi-elastic limit (a — > 1) possibly coupled to large-c limits is certainly of 
interest for inelastic gases 0, 0| . 

The transition rate (|2.2[1 implies that the scattering laws are in all models the same as for inelastic hard spheres, 
i.e. as in an inelastic collision of two perfectly smooth spheres without rotational degrees of freedom, where G and 
g± are conserved, and g'^ = —ag\\ is reflected and reduced in size by a factor a. This implies for the direct collisions, 
(v,w) -> (v*,w*), 

v* = v*(a) =v- i(l + a)p||n 

w* = w*(a) = w + i(l + a)0||n. (2.3) 

The corresponding energy loss in such a collision is, AE = ^(v 2 +w 2 — v* 2 — w* 2 ) = pq\g\\ | 2 , where q = l— p = i(l — a) 
measures the degree of inelasticity. The restituting velocities, describing (v**, w**) — * (v, w), are given by the inverse 
transformation of (|2.3|l . v** = v*(l/a) and w** = w*(l/a). 

A second aspect of the interaction dynamics is the collision rate, K <~ go{g)- In elastic cases the rate is proportional 
to the differential scattering cross-section <r, which depends in general on the relative speed g = |g|, and on the 
scattering angle, which is related to the impact vector n j^]. For elastic hard spheres a(g) = const, and K ~ g, 
and for Maxwell molecules K = const. For repulsive power law potentials, V(r) ~ r _ ", referred to as soft elastic 
spheres, the collision frequency scales as K ~ go~(g) ~ g", and the exponent ;/ is related to the exponent n in the 
power law potential through v = 1 — 2(d — l)/n |47| . For positive n the exponent ^ is restricted to — oo < v < 1. 
As an inelastic generalization of soft elastic spheres we propose systems with a collision frequency that scales as 
K(9-L>9\\) ~ 5^1.911 | CT = ff^^lsil | CT j where the exponent v parameterizes the energy dependence and o the dependence 
on the angle of incidence 9, defined through gu = g • n = cos 9. After inserting K in (|2.1|l two velocity integrations 
can be carried out using the delta functions in (|2.2|) . and the spatially homogeneous Boltzmann equation reduces to 



d t F(v) = I{v\F) = dw 



-K( g± ,^)F(v**)F(w**) - K( g± , 9ll )F(v)F(w) 

a a 



(2.4) 



Here / (• • • ) = (1/0<j) / ^ n (' ■ ■ ) is an angular average over the surface area fid = 2ir d / 2 /T(^d) of a d-dimensional unit 
sphere. We have absorbed constant factors in the time scale. In the one-dimensional case / — > 1 and g — > gn , g± — > 0. 

To clarify the meaning of the exponent a in the collision rate, we consider the impact parameter b, defined as 
b = |g x n| = sin 9 and db — cos 9d9. The distribution V(b) of impact parameters can be obtained from, 



/" dnl^nT ~ / ' d9(sm9) d - 2 \ cos 9\ a - / dbb d - 2 (l - b 2 ) (a - 1)/2 = f dbb d - 2 T(b). (2.5) 
J Jo Jo Jo 

This shows that ~ (1 — fe 2 )( CT_1 )/ 2 for cr < 1 is biased towards grazing collisions, and for a > 1 towards /lead on 
collisions. The cases with cr =^ 1, correspond to pre-collision velocity correlations between the colliding particles. This 
is a violation of Boltzmann's basic postulate of molecular chaos. For dimensions d > 1 molecular chaos requires a = 1 
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or K ~ p^lcos^l, corresponding to a uniform distribution V(b) = 1. For mathematical convenience it looks attractive 
to set a — v 1 and have the simple collision frequency K = \g\\\ v - Then the Boltzmann equation takes the simple form 

M 

d t F(v)=J J dvt\g\\\ v [a- v - 1 F{v**)F{w**)- F{v)F{w)] . (2.6) 

Equation l|2.6(l becomes quite simple for inelastic Maxwell models (v = 0), where all angular dependence in the 
collision frequency is absent [Til \WL l26|. The Boltzmann equation with K ~ |cos#| for the Maxwell model with 
molecular chaos, has been studied in Refs. [sUflL . 

Regarding the subject of interest in this paper, i.e. high energy tails of F(v , t) in spatially homogeneous systems, 
the differences between the classes of models with different values of a are expected to be mostly of a qualitative 
nature |15|. at least for v > 0. As it turns out the shape of the distribution function, in particular its high energy tail, 
depends in a sensitive way fl6L f]~8L IT9 . 26] on the energy dependence of the collision rate at large energies, and not that 
much on the scattering angle. In the present paper the Boltzmann equation 12. 4|) is studied for general exponents 
v and <7, where cases with a ^ 1 do not introduce any additional complications. Most Monte Carlo simulations 
(see Ref. [ill) and most analytic studies [l3L ITil [2(| have been carried out for models with biased distributions, 
V(b) ~ (1 — 6 2 )( <T ~ 1 )/ 2 , and only a few studies exist for molecular chaos model with a = 1 or V(b) = 1 j8lllal48|. 

Equation (|2.4|l is a generalization of the elastic Boltzmann equation to a general class of inelastic models with 
collision rate K ~ g v \ cos0| CT , which include all models presently known in the literature. For v = 1 (n = oo) one has 
inelastic hard spheres and v = corresponds to the softer inelastic Maxwell models (n = 2(d—l)). Negative values of 
v, as n decreases further, are also possible. They would correspond to still softer interactions. Inelastic models with 
v > 1 have also been studied. In the Boltzmann equation with elastic scattering laws there also exists a stochastic 
scattering model, the Very Hard Particle model with v — 2, for which the two-dimensional homogeneous NLBE has 
been solved exactly as an initial value problem |47| . So, there is no a priori mathematical reason to impose restrictions 
on the values of v. Regarding the a exponent we require that the mean collision rate remains bounded. This implies 
that the angular average appearing in the loss term of l|2.4|l should converge. This imposes the restriction a > —1 on 
the models in 1)2. 4|) (see l|A.3[) in Appendix A). Velocities and time have been dimensionalized in terms of the width 
and the mean free time of the initial distribution. Moreover, the Boltzmann equation obeys conservation of mass and 
total momentum, but the average kinetic energy or granular temperature, T ~ (v 2 )t — v o(t), decreases in time on 
account of the dissipative collisions, i.e. / dv(l, v, v 2 )F(v,t) = (1, 0, hdv^t)), where vo(t) is the r.m.s. velocity or 
width of F(v,t). 

To reach a spatially homogeneous steady state, energy has to be supplied homogeneously in space. This will be 
done here by connecting the system to a thermostat or heating device, as discussed in the introduction. For instance, 
a negative friction force may be used as a heating device to compensate for the dissipational losses of energy. Complex 
fluids (e.g. granular matter) subject to such forces can be described -in between collisions- by the microscopic 
equations of motion for the particles, — Vj, and Vj = a.; + £j (i = 1, 2, • • • ), where is a possible conservative or 
friction force per unit mass, and ^ a random force. If the system is driven by Gaussian white noise (WN), the forcing 
can be represented by adding a diffusion term |2§|. —Dd 2 F, to the Boltzmann equation. A negative friction force can 
be included into the Boltzmann equation by adding a force term (d/dv)(a.F). The Boltzmann equation for system 
driven in this manner reads, 

d t F{v)+FF = I{v\F), (2.7) 

where the source term takes the form, 

TF = d ■ (aF) - Dd 2 F = 7o <9 • (vv e F) - Dd 2 F. (2.8) 

Here d = d/dv is a gradient in v— space, and 70 and D are positive constants. When energy is supplied at a constant 
rate, driven dissipative systems can reach a NESS. As expected, the elastic case, driven by white noise, does not reach 
a steady state|23j. 

In the free cooling case the inelastic interactions decrease the kinetic energy of the colliding particles, the r.m.s. 
velocity Vo(t) decreases, and the velocity distribution F(v) — > <5(v) as t — > 00. This scenario is supported by the 
nonlinear Boltzmann equation 12. 411 . because the distribution <5(v) is invariant under the collision dynamics, 

I{v\5) = 0, (2.9) 

as shown in Appendix Al. In general F(v) will be non-Maxwellian, and its shape will depend on the model parameters 
v and cr. The effects of the energy dependence of the collision rate K ~ g v on the high energy tail of the scaling form 
/(c) can be understood intuitively as follows. A tail particle with v 3> vq has a collision rate K ~ |v — w|" ~ v v . 
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The smaller z/, the smaller the value of K in the tail, and the slower the tail particles loose and redistribute their 
energy over the thermal range v < v , and the slower the tail shrinks as a function of v, when compared to the bulk 
velocities. This leads to an increase in over-population of the tail. The reverse scenario applies when increasing v, 
leading to a decrease of over-population. 

Intuitive pictures of the effects of the heating devices can also be developed. The linear Gaussian thermostat with 
a = 70V simply produces a linear rescaling of the velocities, and has the nice property of supplying energy to or 
subtracting energy from the system without changing the scaling form /(c). The nonlinear friction force a = 70 ~vv e 
produces a nonlinear rescaling, A(v) v. For 6 > 1 it represents a heating device that puts selectively energy into the 
tail particles, thus leading to an increase of over-population of the tail. For 6 < 1 the reverse scenario applies. The 
white noise forcing on the other hand adds randomly velocity increments Aw, drawn from a uniform d— dimensional 
distribution, to thermal and tail particles. So, it is more efficient in redistributing and randomizing velocities of tail 
particles, especially at smaller v— values. 

If the collision frequency K oc g" gets smaller (c.q. larger) at large relative velocities, then the tail distribution 
shrinks slower (c.q. faster) than a Maxwellian, resulting in an overpopulated (c.q. underpopulated) high energy tail, 
described by a stretched (c.q. compressed) Gaussian, F(v) ~ exp[— [3v h ] with < b < 2 (c.q. b > 2). The limit 
as b — > + corresponds to heavily overpopulated power law tails. Similar scenario's presumably occur also when 
dissipative systems are coupled to an energy source, and are slowly heating up or reaching a non-equilibrium steady 
state (NESS). 

To study the behavior of the v.d.f. with a shrinking width we consider the distribution function /(c), rescaled by 
the r.m.s. velocity vq, i.e. F(v,t) = Vq f(v/vo). Suppose that for large c = v/vq the velocity distribution can be 
separated into two parts, /(c) = /o(c) + h(c), where h(c) is the singular tail part, and fo(c) the presumably regular 
bulk part 49]. The tail part h(c) may be exponentially bounded ~ exp[— (3c b ] with < b < 2, or of power law type. 
In the bulk part /o(c) the variable c is effectively restricted to bulk values in the thermal range v ^ vq or c < 1. 

To obtain an asymptotic expansion of /(c) we consider the ansatz, 

f(c) = S(c) + h(c), (2.10) 

and we linearize the collision term I(c\f) around the delta function. This defines the linearized Boltzmann collision 
operator A as, 

I{c\6 + h) = -Ah(c) +0(h 2 ), (2.11) 

where Ea. l|2.9|l has been used. Here we restrict ourselves to isotropic functions h(c), depending only on c = |c|. As 
will be shown in the next section the eigenfunctions of A decay like powers c~ s . Consequently they are very suitable 
for describing power law tails, /(c) ~ c~ s . 

In case the high energy behavior of the v.d.f. 's are stretched Gaussians ~ exp[— /3c b ], the above eigenfunctions are 
no longer useful, and we have developed a method to determine an asymptotic expansion of the form, 

In /(c) - -f3c b + f3'c b ' +xlnc+-- - , (2.12) 

as will be discussed later. 

A comment seems in order here. The delta function in l|2.11|l does not mean that the thermal bulk part /o (c) ( with 
c > 1) is singular, as explained in the paragraph above l|2.11|l . In fact, it is irrelevant for our arguments of extracting the 
dominant large-c singularity whether the thermal part may or may not contain any singularities. The separation of f(c) 
into a 'regular' part /o(c) and a singular part h{c) is completely analogous to the method used in Refs. mill in m 
for predicting the power law tail for Maxwell models using the Fourier transform method. Let the Fourier transform 

of /(|c|) be T k f{\c\) = f(k). Its small-* behavior takes the form f(k) = 1 + k 2 f 2 + k 4 f 4 H h Ak a H with 

a 7^ even integer. Here the series expansion in powers of k 2 = \k\ 2 represents the regular part, Tkfo(c) = fo(k) 
1 + fc 2 / 2 + fc 4 / 4 + • • • , whereas Ak a represents the dominant the small fc 2 -singularity of h(k). Subsequent Fourier 
inversion of the small-fc behavior of f{k) shows that the regular bulk part /o(c) can be viewed in zeroth order of 
approximation as 5(c). For c > 1 only the tail contribution of h(c) = /(c) — 6(c) survives. The dominant small-fc 2 
singularity of h(k) corresponds to the dominant large-c singularity of h(c). 

III. SPECTRAL PROPERTIES OF COLLISION OPERATORS 

The first part of this section is devoted to a study of the eigenvalue problem of the linearized Boltzmann collision 
operator A. This is most conveniently done by first constructing the adjoint operator, A^, and determining its 
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eigenvalues and eigenfunctions, also referred to as (left) L-eigenfunctions of A. The adjoint is defined through 

(k\Ah) = (h\A*k), (3.1) 

where the inner product is (k\h) — J dck(c)h(c), and k(c) and h(c) are isotropic. This can be done by using the 
relation, 

(k\I(F)) = \J J dcdwK(g ± , gil )F(c)F(w) [k(c*) + k(w*) - k(c) - k(w)} (3.2) 

with g = c — w, as follows from 12.4fl . To facilitate the analysis we express c*, w* in Eq. (|2.3[1 in terms of cj_, wj_, 
cm, wu with coefficients p = |(1 + a) and q = |(1 — a) (see Appendix Al). Substitution of l|2.10() into (|3.2() . and 
linearization yields 

(k\I{6 + h)) ~ -(k\Ah) = -(h\A^k) 
= - /„ / dcK{c±,c\\)h{c) [k(c) + k(0) - k(\c± + nqc\\ |) - fc(|npc„ |)] , (3.3) 

where S(w) has been integrated out. As fc(c) = fc(|c|) is isotropic, the adjoint becomes 



A f fc(c) = / K(c ± , C]] ) [fc(c)+fc(0)-fe(c|c x +ngc|||)-fc(pc|c|||)] (3.4) 

J n 

with K (cj_, cii) = c^ _(T |c|| |°". In one dimension, where c_l = and f — > 1, the adjoint takes the simple form 

Atfc( c ) = \c\» [k(\c\) + k(0) - k(q\c\) - k(p\c\)} . (3.5) 

Inspection shows that insertion of k(c) = \c\ B generates a new power of |c|, which solves our eigenvalue problem 

A^\c\ s = \c\ s+ "[l + S s , -q s -p s ] , (3.6) 

where A s = 1 + 5 8i o — q s — p s is the eigenvalue, valid for s > 0. Strictly speaking, |c| s is not an eigenfunction of 
but Eq. (|3.6|l shows that the family of power law functions is stable under the action of At. In the following, we 
will employ the terminology "eigenvalue" and "eigenfunction" with a similar slight abuse of vocabulary. The same 
argument shows with the help of the homogeneity relation K(c±, cy) = c v 'K(c±, C||), that k{c) = c s with s > is also 
an eigenfunction of A^ for dimensions d > 1. The eigenvalue follows as 



A s (<7 ) = / | C || | ff [1 + S a<0 - \c± + nqc {] \ s - p s \c {l \ s ] . (3.7) 

■/ n 

For fixed c the integrand of this (d — 1) dimensional integral is a function of cy = c ■ n = cos 8. It reduces to a single 
integral over a polar angle (0 < 8 < tt) for d > 3, and for d = 2 to an integral over an azimuthal angle (0 < 8 < 2tt). 
The integral is evaluated in Appendix A2, with the result valid for s > 0, 

X s (a)=f3 (T {l + S sfi - 2 F 1 (-f,£±i;2±£i|l-g2)} _p^ a+(T . (3.8) 

Here 2^1 (a, 6; c|z) is a hyper-geometric function and j3 s = f \a ■ n| s is an average of | cos^| s over a solid angle, given 
in i|A.3|) . It is well defined for a > — 1. 

So far we have constructed in (|3.6|l - (|3.7|l the eigenvalues A s (<r) and corresponding eigenfunctions k s (c) = c s of the 
adjoint At, which are the (left) L-eigenfunctions of A. To determine the (right) R-eigenfunctions h s (c), corresponding 
to A s (a) , we need to construct A in (|2.11l) explicitly The steps are rather technical, and are carried out in Appendix 
A3. However, the procedure is quite similar to the steps from (|3.4(l to Ij3.7|l . Once A is constructed, inspection shows 
again that Ac 1 * = /i r c r+l/ , where the new eigenvalue fi r is different from X s . Next, r in /i r is chosen such that /j. r = X s . 
This yields r = —s — d — v. 

The most important spectral properties for our purpose are: 
(i) For s > the eigenvalues and R- and L- eigenfunctions are ^l] (with details derived in (|A.16|I and l|3.6M3.7p ). 

Ac- s - d -» = X s (a)c- s - d 

AV = X s (<j)c s+v . (3.9) 

For s = one has the stationary eigenfunctions, which are invariant under collisions, 

A<5(c) = and A f ■ 1 = (A (cr) = 0) . (3.10) 



7 



R- and L-eigenfunctions are different because A is not self-adjoint. There is in fact among the eigenfunctions in i|3.9|) 
another, less trivial, stationary R-eigenfunction, c~ a ~ d ~ v , where s = a* is the root of X s (<r) = (see Fig l3.1|) . 
(ii) The spectrum exists for all s > 0. It is continuous for s > 0, and s = is an isolated point of the spectrum. At 
s = one has 2-Fi(0, b; c\z) = 1 and consequently 



Ao(ct) = 



and 



lim A s (cr) = -/3 CT 



(3.11) 



The zero eigenvalue expresses the conservation of mass in binary collisions. In the elastic limit A2(er, a = 1) = 0, 
which expresses the conservation of energy (see Fig l3.1|) . 

(hi) As Ps+a, given in (|A.3(1 . and Ai (<r), given in (|A.4(I . are monotonically decreasing with s, A s (er) is monotonically 
increasing with s, and approaches 



lim A s (cr) = (3cr. 



(3.12) 



(iv) A s (er) is independent of v, i.e. the energy-dependence of the collision rate K ~ <7 1 '|<?||| CT does not affect the 
eigenvalue. So, it is the same for inelastic Maxwell molecules, hard spheres, very hard particles, and very weakly 
interacting particles, as modeled by Eqs. 1|2.1|) - i|2.4[l . The reason is presumably that the scattering laws are the same 
in all models, and equal to those of inelastic hard spheres (|2.3|) . 

(v) The eigenvalue \ s {o~) does depend on the angular exponent a. We also recall that A s (cr) for the molecular chaos 
models has a = 1. For the case (d = 3, a = 1) the eigenvalue reduces to a simple expression, 



A s (a = l,d = 3)//3 1 = 1- 



s + 2 



(!-<?)* + 



1-9 



s+2 



1 



(3.13) 



obtained before in Refs. [3,[l5| for the three-dimensional Maxwell model (y = 0, a = 1), satisfying molecular chaos. 
To obtain this result form (|3.7(l one may use the relation, 



1-(1 



Z) 1+ 1 S 



(3.14) 



One should also keep in mind that the ^-dependence of A s (a = v) in the simple model of Eq. (|2.6|) only refers to the 
angular dependence of the collision rate K = \g\\\ v — g v \g\\\ u . 

(vi) The eigenvalue for s = 2 has for all d the simple form 

A 2 (a) = 2 Pq p a+2 , (3.15) 
as follows from the relation 2-F\(— 1, b; c\z) = 1 — bz/c and (|A.3|I . 

(vii) The behavior of \ s (a)/(3 a as a function of s is shown in Fig l3.1l for various values of the inelasticity, where a = 1 
and d = 2. 

Before concluding this section we analyze the asymptotic form of the collision operator I(c\f) at large c, in order 
to improve the estimate, Ioa(c\f) ~ —f3 a c u f, given in the literature |18| . This extension is needed to obtain an 
asymptotic expansion of the form (|2.12|) . or equivalently, 



/( c )~c*exp[-/3c 6 +/3'c b '], 
with b > b'. The analysis is given in (|B.9|) of Appendix B, with the result for c ^> 1, 

hb(a + l) 



-/oo(c|/) = Aoo/(C) ^ PaC» 



l-JC a c~T 



(3.16) 



(3.17) 



The coefficient is given by 



/C 



3£) / 2 \ 
SET) [(iba-q 2 ) J 



£1 

r( 

(d-l) 
/3fc(l-9 2 ) 



(a+l)/2 



(d^l; general) 

(cr = l; molecular chaos) 



(3.18) 



These relations are valid for all x,/3', b' < b, and d > 1. For c? = 1 there are only exponentially decaying corrections, 
and JC a = 0. The arguments, used here and in the appendices, to analyze the properties of the operators A and 
are similar in spirit to those used in Ref. |32| . 
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IV. SCALING EQUATIONS 

As sketched in a qualitative manner in Section II, the velocity distribution function F(v, t) shrinks -without external 
driving and in zeroth order approximation- to a delta distribution <5(v) as the thermal velocity decreases. However, 
if one re-scales the velocities in the shrinking distribution in terms of its instantaneous width Vo(t), then the rescaled 
distribution rapidly approaches a time-independent scaling form /(c), 

F(v,t) = (v (t))- d f{c) and c = v/v Q (t), (4.1) 

as discussed in Refs. jHIISSil. Direct Monte Carlo Simulations (DSMC) of the Boltzmann equation [T l ITfl l El 
E3 | have confirmed that after a short transient time the rescaled velocity distribution can be collapsed on a time- 
independent scaling solution /(c). These observations indicate that the long time behavior of F(v,t) in freely cooling 
systems approaches a simple, and to some extent universal, scaling form /(c), which is the same for a general class of 
initial distributions. It satisfies the normalizations, 

J dc{l,c 2 }f(c) = {l,±d}, (4.2) 

consistent with the definition of J dvv 2 F(v,t) = \dv^(t) of the r.m.s. velocity Vo(t). 

To study the scaling forms in the case of white noise driving, we consider the equation for the mean square velocity, 
using the Boltzmann equation (|2.7[) - (|2.8|) . 

d A = l div 2 lI{F))+ 2 D{v 2 ld 2 F) 

= -27(1/, cr)<+ 2 + AD. (4.3) 

Here, we have used ()3.2|) with k = v 2 , and the relation AE = —pqg 2 below Ij2.3|l . performed the rescaling in (|4.1|l . 
and introduced 

7(1/, a) = i TO /W(|c- Ci r 2 )). (4.4) 

The average ((...)) is calculated with the time-independent weights, /(c)/(ci). So 7(1/, a) is an unknown constant, 
except for the cases v = (Maxwell) and v = — 2 (WN-threshold model), where it takes the values, 

7(0,(r) = pq {3 a+2 and a) = pq -j f3 a+2 . (4.5) 



For the case of free cooling without energy supply (D = 0), the r.m.s velocity in Ij4.3|l keeps decreasing for large t, as 

v (t) 



«o(0) 



(l + vt/t^v,*))- 1 /" (z/>0) 

exp[-t/* c (0,<7)] (i/ = 0) , (4.6) 

{l-\v\t/t c {u,a)Y'\"\ (v<0) 



where the constant t c (v, o~) = <j)vq(0)] is unknown except for v = {0, —2}: 

t c (0,a) = l/pqpa+2 and t B (-2,a) = d/pqp a+2 . (4.7) 

Equation 14.6|l is the extension of Haff 's homogeneous cooling law [50| to inelastic soft spheres. 

The most important time scale in kinetic theory is the inverse of the mean collision rate, defined as the average of 
the collision rate K = g v \g\' J : i.e. 



w(t) = / / dxdw g v \g\\\ a F{v,t)F{w,t) 

= l{{\ Cl -cn)v v Q {t)=tM<{t)- (4-8) 

The second line above applies only to scaling solutions, where Q{v, a) is known explicitly for C(2, a) = d[3 a (Very Hard 
Particles) and C,{Q,a) — j3 a (Maxwell). The collision frequency u>(t) depends on the external or laboratory time t. 
We also consider the collision counter or internal time of a particle r, defined through the relation dr = u>{t)dt. It 
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represents the total number of collisions r that a particle has suffered in the external time t. With the help of i)4.8)) 
and 1)4.6)) this differential equation can be solved to yield 



1 



a) 



In 



l + v 



t c (v, a) 



l/u 



where 



-cu(v, a) 



7(1/, a) _ TO /3 ff+2 «|c-c 1 r+ 2 )) 



CM 



^((Ic-dh) 



is an unknown constant, as it depends on /(c). It reduces for v — to tn(0, a) = pq/3 CT +2/ 'At- By combining 
1)4.6)1 in the case of free cooling, the homogeneous cooling law takes the simple exponential form 

Vo(t) — vq(0) exp[— w(v, <j)t\. 



(4.9) 

(4.10) 
and 

(4.11) 



On the other hand for v < the r.m.s. velocity vo(t) in (|4.f>|) . and the granular temperature v^it) vanish in a finite 
time t c (v,<j) /\v\ — 1/7(1/, ct)uq(O). At the same time the collision counter r in 1)4. 9fl diverges. This means that in 
inelastic soft sphere models with v < an infinite number of collisions occurs in a finite time t c (i>, a) /j v\ . A collision 
counter, diverging within a finite time, is also observed in the phenomenon of inelastic collapse 51], occurring in 
inelastic hard sphere fluids. There a finite number of particles line up in a linear array in configuration space. The 
divergence of the collision counter r in (|4.9|l at t — t c (y^ a)j\v\ is a reflection in velocity space of this inelastic collapse 
phenomenon. In the WN-driven case one may also analyze the collision counter r through dr — u(i)<it, and one 
observes similar phenomena of inelastic collapse for v < — 2. 

Next we consider the NESS. In the case of WN-driving the energy balance in 1)4.3)1 has a stationary or fixed point 
solution ?;o(oo), given by 



u 2b {oo) = 2D/j(v,a) (2b = v + 2), 



and the energy balance equation 1)4. ?>\ can be expressed as, 

dvl{t) 



dt 



AD 



vp(t) 

V (oo) 



2b' 



(4.12) 



(4.13) 



The solution vo(t) may approach to or move away from vq(oo), depending on the sign of b. For 6 = l + ^/2>0, 
the fixed point solution is stable and attracting, for b < it is unstable and repelling, and for b = it is marginally 
stable. In a system at the stability threshold (6 = 0), a marginally stable NESS only exists if one can fine tune the 
driving parameters in (|4.3|l to D — j(— 2, u)/2 — pq(3 a+ 2/d- Otherwise, the energy for 6 = or v — —2 is increasing 
or decreasing linearly with time, 



vl{t) = «g(0) + (AD - 2 P q/3 a+2 /d)t. 



(4-14) 



The possible existence of time-independent scaling solutions can be investigated by substituting 14.1fl into l|2.4|l with 
the result 



D 



d 2 f 



7(1/, a) 



2D 

„,i/+2 



D 



(4.15) 



On the last line vq has been eliminated using (|4.3|) . This equation is not yet an integral equation for a scaling function, 
because the coefficients still depend on time through vo(t). 

There are two possibilities for time-independent solutions. The first one is obtained by setting D = 0. This gives 
the integral equation for the scaling solution in the free cooling state, which has been studied extensively in the 
literature |l6|, Il8l |26j , and leads to high energy tails of stretched exponential form for v > or to power law tails for 
v = 0. The second possibility for a time-independent solution is for the stable NESS in the WN-noise driven case 
with i>o(oo) given by 1)4.12(1 . Then the coefficient of the first term on the r.h.s. of 1)4.15)1 vanishes, and the integral 
equation reduces to, 



I(c\f) = -i 7 ( V) a) d*f = -1 7(1/, a) (/" + ±±f) 



(4.16) 
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where a prime denotes a c— derivative. 

However, at the stability threshold (y = —2) there exists a marginally stable state, described by l|4.14|) with 
vq(oo) = v (0), and obtained by fine tuning the driving parameter to the value D = \"f{— 2, a). This yields the scaling 
equation for the WN-noise threshold model with v = — 2, 

Kc\f) = -^pq^ +2 d 2 f = -1a 2 ( ( t)9 2 /, (4.17) 

where <|4.5|) and i|3.15[l have been used. Furthermore, it is interesting to note that the existence of a time independent 
scaling form does not require the energy v 2 to be stationary in the WN-threshold model (y = —2). Suppose we add 
an additional friction force a = 7ov, the Gaussian thermostat, to the source term of the Boltzmann equation (|2.7|l . 
(where 70 may be positive or negative), then the energy balance equation (|4.3|) and the integral equation for the 
scaling form 14.15fl both get an extra term, i.e. 

voVo = -7(- 2 ) o-) + 2D + 7 Uq 
I(c\f) = [-v vo+lov 2 }d-(cf)-Dd 2 f 

= [ 1 {-2,a)-2D]d-(cf)-Dd 2 f. (4.18) 

When vq is eliminated from the first two equations, both extra terms containing joVq cancel, and the resulting equation 
is identical to H4.15[) with v = — 2. This is an integral equation for the scaling form. 

Subsequent fine tuning of the diffusion coefficient to the value D = \^{— 2, a) simplifies the integral equation to 
the same scaling equation (|4.17|l . and the energy balance equation to vq — 70 «o- Consequently in the limit of large 
time energy may diverge or vanish, rather than stay constant. So, the total energy does not have to be finite for the 
existence of a scaling form, and the time independent scaling solution /(c) of the fine-tuned WN-driven Boltzmann 
equation in (|4.17|) is not affected by adding an extra Gaussian thermostat. 

In performing DSMC simulations it should be noted that the realization of a marginally stable NESS is more 
complicated for driving by a stochastic force than by a deterministic friction force. The reason is that D is the mean 
strength of the random kicks (£ 2 ), and the realized strength fluctuates around the mean. So fine tuning the value of 
D by choosing it to satisfy (|4.12l) . does not make vq vanish exactly, but it may be a bit positive or negative. Then by 
re-scaling the velocities respectively down or up at regular intervals, i.e. by using the Gaussian thermostat, a = —70V, 
respectively with 70 > one can decrease v to 0, or with 70 < increase v to 0, thus reaching the NESS. This 
addition of an extra Gaussian thermostat leaves the integral equation for the scaling function invariant, as we have 
seen in (|4.18|l . One may also consider an alternative but very similar algorithm, as given in Ref. |25| . 

Exact closed form solutions of the scaling equations are not known, except for the freely cooling one-dimensional 
Maxwell model ^l]]. However, the high energy tails of /(c) can be calculated explicitly with the new methods, 
presented in this article, both for exponentially bounded tails, exp[— /3c b ] appearing in stable NESS, as well as for 
power law tails, /(c) ~ c~ a appearing in marginally stable NESS (b = 0). 



V. ASYMPTOTICS FOR STABLE NESS (y > 2) 

In this section we focus on the high energy behavior of the scaling form /(c) for a stable NESS (6=1 + v/2 > 0) 
for inelastic soft sphere models as determined by the integral equation (|4.16l) . For large c- values the collision operator 
reduces to I{c\f) ~ — Pg c" f (c) and the solution of (|4.16(l has the form /(c) ~ exp(— (3c b ) with 6=1 + v/2, as has 
been derived in Refs fl8L l26j | . The goal of this section is to determine the sub-leading correction in the asymptotic 
expansion (|5.1|) of In /(c) for large c, i.e. 

ln/(c) ~ -(3c b +(3'c b ' +xlogc + logA + ... (5.1) 

where 6 > 6'. To determine these exponents and coefficients we need a better asymptotic estimate of the Boltzmann 
collision operator 7(c|/), when acting on a function of the form (|3.16(l . 

According to Appendix B the large— c behavior of I(c\f) is determined by the linear operator Aoo, which is given 
in IjB.lfl . It is a limiting form for c 3> 1 of the linearized Boltzmann collision operator A, and its action on functions 
of the form (|3.1t)|) has been analyzed in Appendix B with the result l|B.9|) . 

Ioc(c\f) = -Aoo / ~ (3 a c v [1 - K a c~ a ] /(c), (5.2) 

where JC a is a constant given in l|B.8() and a = 6(er + l)/2. Then the scaling equation (|4.16l) becomes 



fl, c" [1 - K a c- a ] f = §7(1/, a) {/" + (**=i) /'} 



(5.3) 
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As this equation is homogeneous, the constant A in l|5.1|l cannot be determined. So we set A = 1 and for consistency 
impose the restriction b > 6' > 0, as c~' b ' <C | In for c ^> 1. The parameters in (|5.1|1 can be determined by 
substituting /(c) of (|3.16|) into 1|5.2[1 and matching the exponents of powers of c and equate the coefficients, i.e. 

r) a c v [1 - K a c- a ] ~ /3 2 6 2 c 2b - 2 - 2/3bl3'b' c b+b '- 2 - f3bc b ~ 2 [d - 2 + b + 2 X ] + ... (5.4) 

To leading order we equate the terms proportional to c v and c 2b ~ 2 , which yields the exponent b and the coefficient (3, 
with the well known results 261, 



/(c) ~ exp[-/?c b ] , 6=l + ±^>0 , bp = Va = ^2jj^a). (5.5) 

The dominant tails are stretched Gaussians (0 < b < 2), or compressed ones (b > 2). These tails are respectively over- 
populated or under-populated when compared to a Gaussian. The results for soft sphere models with b = 1 + \v > 
are valid for all dimensions d > 1. We also note that the exponents 6 found here are equal to the b- values that 
determines the stability threshold in the energy balance equation (|4.13|) . In principle, <x) in (|5.5|) can be calculated 
perturbatively, using the method developed in Ref [18, 26], or it can be measured independently using DSMC, as will 
be done in the present paper. 

We first consider Eq. (|5.4|l for the molecular chaos models with a = 1, where (|lj.8(l gives (3bK\ — (d — 1)/(1 — q 2 ), 
and a = b = 1 + \v. The only consistent matching of the sub-leading exponents (satisfying the restriction b' > 0) is 
obtained by setting j3' = 0. The parameters in the sub-leading correction are then 

/3' = , X = -\{d-l)-\v+\m 1 = y0£-\ V (5.6) 

valid for all allowed q- values Q < q — \ (1 — a) < i. For a ^ 1 one can easily verify that the exponents in l|5.4|l obey 
the relation v — a ^ b — 2, and matching gives either v — b = b + b' — 2 > b — 2 or the reversed inequality. In the 
former case, b' = 6(1 — <j)/2 > which is realized for a < 1. Then equating coefficients yields (3' = piC a /(l — a). In 
this case the sub- leading asymptotic correction in (|5.1|) is (3'c b \x\ logc, and in the spirit of asymptotic expansions 
we set x — 0- 

For the reversed inequality b' = 6(1 — <j)/2 < 0, realized for a > 1, the sub-leading asymptotic correction in l|5.1|) 
is |x| logc ^ (3'c b , and in the same spirit as above, we set (3' = 0. The value of exponent x 1S found by setting the 
coefficient of c b ~ 2 in l|5.4|l equal to zero. In summary: 

£'=0, X = -\{d- 1 )-\v ( CT>1 ) ( 5 - 7 ) 

b' = \b{l-a), F=pK a /(l-a), X = 0, (\a\ < 1) (5.8) 

Next we focus for the WN-driven case on comparison of the Monte Carlo simulations with the theoretical predictions. 
Fig l5. II shows the DSMC data of one-dimensional soft sphere models for /(c) as a function of c b with b = 1 + for 
several values of v. Moreover, the asymptotic large— c prediction, ln/(c) ~ — (3c b ~ (v/i) lnc in l|5.6l) . is confirmed in 
Fig l5.2l The value of (3 is independently obtained by measuring 7(1/, a) in (|5.5(l during the DSMC simulation. Note 
that for the one-dimensional Maxwell model the prediction /(c) ~ exp[— cy / 2/pq] has a very simple form, because 
(|5.5|) gives 770 = y^2/pq as (3 S = 1 for all s. Moreover, this figure shows that the correction, — (^/4)lnc, leads to an 
improved agreement between theory and simulations, even though it is tiny because of the limited range of velocities 
accessible. 

The two-dimensional WN-driven model, with a = v = —1/2, is illustrated in the next two figures. Fig l5.3l first 
displays the measured scaling form /(c), plotted versus c b with 6 = 1 + v/2 = 3/4 for different values of a. The straight 
parts of the high energy tails confirm the dominant velocity dependence, In /(c) ~ — /3c 3 / 4 , with (3 increasing with a. 
For large enough a, /(c) is getting closer to a Gaussian. Moreover, Fig l5.4l presents a very interesting illustration of 
the sub-leading correction exp(/3'c b ). While a plot of In /(c) versus c b yields a straight line, and thus an apparent 
agreement with the dominant theoretical prediction exp(— {3c b ) : closer examination shows that measuring (3 through 
a fit to exp(— (3c b ) would disagree with the predicted value of (3 in (|5.5|l . which is computed from the numerical data 
using the definition 15.5fl which involves 7(1^, a). On the other hand, comparison of /(c) with the full expression 
exp(— j3c b + (3'c b ) -where (3' is measured from its definition (|5.8Jl involving JC a - gives a very good agreement. 

Regarding the soft sphere systems in stable non-equilibrium steady states with model parameters v > — 2, we may 
conclude that the agreement between analytic and DSMC results for high energy tails is very good. Essentially all 
analytic and numerical results in this section are new. They are the leading and sub-leading asymptotic correction 
factors in the scaling function /(c) in a stable NESS (6 > 0) for the general class of inelastic soft sphere models. These 
results include the few special cases known in the literature, i.e. the dominant asymptotic form /(c) ~ exp[— [3c b ] 
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for white noise driven inelastic hard spheres (y = 1) and similar results for inelastic Maxwell models (y = 0). 
The only known result about sub-leading correction factors concerns the very special case of the three-dimensional 
molecular chaos Maxwell model (d = 3, v = 0, a = 1), for which Bobylev and Cercignani have derived the result 
/(c) ~ c x exp[— (3c b \ with \ = <? 2 /(l — Q 2 )- Although the derivation is not very transparent, their result agrees with 
(|5.6() for this special case. 



VI. POWER LAW TAILS AT MARGINAL STABILITY {u = 2) 

Marginal stability is a limiting property of a stable NESS as b — > + , which can not occur in free cooling, but only 
in driven states, here driven by white noise. Below l|4.14|l we have explained that the marginally stable NESS can not 
be reached by natural time evolution, but only by fine tuning the external parameters in the energy source term. 

As we have seen in l|5.5fl the asymptotic solution /(c) for stable states (b = 1 + \v > 0) is to leading order described 
by /(c) ~ v4exp[— /3c b ] with = r\ a jb. To illustrate how power law tails arise, it is convenient to set A — e@, and take 
the limit of /(c) as b — > 0. The result is /(c) ~ exp[— n a (c b — l)/b] — + c~ n ° , with rj a = limf,_o Va- Of course r\ a is not 
the full exponent of the tail, because the exponential form above represents only the leading asymptotic behavior for 
b > 0. For instance, any correction factor exp[— (3'c h ] as appearing in (|3.16|1 . where b' = B(b) — > as b — > 0, would 
give additional contributions. 

To determine the exact exponent of the power law tail we have to solve the scaling equation (|4.17() for the WN- 
thrcshold model, where v = —2. As we expect power law solutions from the arguments above, we use the linearization 
in JZHl-lEinj to give, 

-I(c\5 + h)~Ah = ±-\ 2 (a)d 2 h = ±- pq f3 a+2 d 2 h, (6.1) 
4a 2d 

where /(c) ~ h(c) for c > 1. Inspection of this equation shows that the operators on both sides of the equation have 
the same set of R-eigenfunctions, c~ s ~ d+2 . Substituting this function into (|6.1|l gives the transcendental equation, 

X s (a) = (l/4d)s(s + d- 2)A 2 (a) = (l/2o>(s + d - 2)pq/3 (T+2 . (6.2) 

The possible roots of this equation, s = a^ct), depend on the spatial dimensionality d, and on the coefficient of 
restitution a through p = 1 — q = i(l — a). They are candidates for power law exponents in an asymptotic solution 
of the form, 

/(c) ~ c- a - d + 2 (6.3) 

for the WN-threshold model. However there exist a priori restrictions on the exponent a through the normalizations 
(|4.2|) . If one is interested in solutions with bounded energy, then the constraint, J dec 2 /(c) < oo, or equivalcntly 
a > 4, applies for all d. If one allows the energy to keep increasing - as discussed below (|4.18|l - one still has to 
obey the restriction, J dcf(c) < oo, implying a > 2 for all d. Moreover, for d > 1 the restriction a > — 1 applies, as 
discussed below (|2.6H . At d = 1 the angular exponent a is absent. 

To start we consider the one-dimensional case, where Ea. H6.2(l reduces to 1 — p s — q s = ^s(s — l)pq. It has two 
solutions: sl = 1 and sr = a — 3. Only the largest root satisfies the normalization constraint, and corresponds to 
the solution, /(c) ~ l/c a+d_2 ~ 1/c 2 . We also note that in the one-dimensional case the power law exponent a is 
independent of the coefficient of restitution. In our DSMC calculations the energy is kept bounded by applying the 
Gaussian thermostat, which does not affect the integral equation for the scaling function, as discussed at the end of 
Section IV. 

This analytic result is in good agreement with the one-dimensional DSMC results as shown in Fig l6.ll where the 
algebraic tail is observed over more than 3 decades in /(c). The energy in this state is however infinite. Why this 
high overpopulation of the power law tail? The WN-threshold model is very inefficient in equilibrating its high energy 
particles in comparison with hard spheres (y = 1), and even with Maxwell models (y = 0), because the tail particles 
rarely suffer a collision as their collision rate decreases as K ~ \c\ v ~ l/|c| 2 . 

To obtain a qualitative understanding of the a— and d— dependence of the exponent ad{a), we use a graphical 
solution method by plotting the l.h.s. y s and the r.h.s. y s in Ij6.2|l as functions of s, and determine the two intersection 
points. Here we have defined, 

Vs = — and ys = u^Td) s{s + d - 2) > (6 - 4) 

where the relation Pcr+2/Pa = (cr + I) / (<r + d) has been used. The procedure is illustrated in Fig 16. 21 for (d = 2, a = 1). 
In the elastic limit (a — > 1~) the pre-factor in y s becomes smaller, and the right intersection point, sr = ad(ct), moves 
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to infinity, ft satisfies the condition a > 4, and is an acceptable power law exponent. The left intersection point, sl, 
remains less than 2, and does not satisfy the constraints of finite energy and mass. 

The numerical solution of (|6.2|1 for (d = 2, a = 1) yields ad(a), shown in Fig l6.3l Fig l6.4l shows the good agreement 
between the predicted power law tail and the DSMC measurements of /(c) for d = 2 and various values of a. 

Guided by the numerical result, 02(0; = 0) ~ 4 , we have verified that this value is an exact solution of 16. 2f) for 
d= 2. To do so we have calculated A 4 (l) using the relation 2 i<i(-2, 1, \{d+ l)|3/4) = (d 2 + d-3/2)/(d + l)(d + 2), 
as can be obtained from the second degree polynomial 2-Fi( — 2, b; c\z) in z, given through the series in i|A.6|l that 
terminates after the term with n = 2. For higher dimensions no such simple solutions seem to exist. For instance, 
ad(a = 0) = 4.929, 5.812, 6.672... for d = 3, 4, 5, .... As adia) is an increasing function of a, our graphical and numerical 
analysis shows that the power law tails, found in the present section for d > 2, all carry a finite energy. 

For general d the transcendental equation l|6.2|l can only be solved in a few limiting cases. By taking the elastic 
limit (q — > 0) in l|6.2|l . and comparing right and left hand sides one sees that the solution s must become large as 
q — > 0. So we need A s (<r = 1) for large s, as given in (|3.12l) . and the dominant behavior of a at small q is, 

= y/2d(d + a)/q(l + (r) {q^0;d>l). (6.5) 

further sub-leading corrections to (|6.5|l can also be calculated using asymptotic results for 2^1(0, b; c\z) at large a. 

For large dimensionality we take the coupled limit as d — > 00 and s — > 00 while keeping x = s/d = fixed, and we 
set cr=l. This leads to an algebraic equation of third order with three roots {x- < 0, xq — 0,x + > 0}, where the 
largest one determines the exponent in the power law tail, and reads, 

a d (a) =x+d~ 2g(1 l g2) {VM(l + q)(l-q 2 ) + q 6 -2q + q 3 } (d-foo). (6.6) 

The behavior of x+ vs a is illustrated in Fig l6.3l The results (|6.5(l and i|6.6[l agree to dominant order in the respective 
limits d — > 00 and q — > 0. 

Finally, we discuss the power law tail generated by an energy source "at infinity" , as recently proposed in Ref. • 
Here the heating device injects energy in the ultra high energy tail of the v.d.f. It works as follows. By randomly 
selecting a particle at a very small rate 7 <C loq (being the mean collision rate), and giving the selected particle a 
(macroscopic) amount of energy AE(j), equal to the total energy lost by the system in the typical (long) interval 
I/7 between two injections. After a short transient time, the energy cascades down to lower and lower energy, and in 
that manner a NESS-v.d.f /(c) is maintained that does not evolve under inelastic collision dynamics. 

The mathematical framework developed in the present paper, is very suitable to discuss the problem posed above. 
We are in fact looking for a solution of equation I(c\f) = 0. This implies through Eqs. H2.10fl - H2.11D that we 
want a solution of Ah(c) = 0. The answer is given by (|3.1U|) . The required scaling function is the i?-eigenfunction 
h(c) ~ c~ s ~ d ~", where s is the root of A s (ct) = with A s given in 1)3. 8f) . This expression is identical to Eq. (11) of 
Ref. [3~2, as can be verified from l|A.ll|) . 

We recall from property (iii) below l|3.12|l that X s (a) is independent of ^-exponent in the collision rate K ~ 
g v \ cos^l 17 , and the er-exponent refers only to the angular dependence. For molecular chaos models a = 1, as shown 
in H2.5I . and we restrict ourselves to this case. 

The graphical solution to this problem can be seen in Fig l3.ll as the point of intersection of the A s -curve and the 
s-axis. For d = 1, the eigenvalue reduces according to l|3.6|l to 1 — p s — q s = with solutions s = a* = 1, and the 
scaling form /(c) ~ c —a—*-v ^ c -2-v ^ -p Qr dimensions d> 2 the scaling form has the power law tail, 

/(c) ~ (6-7) 

where s = a* = a* d {a) is the solution of X s (a = 1) = and depends on d and a. 
From general considerations one obtains the bounds [32 

1 = at(a) < a* d (a) < a* d (l)=2, (6.8) 

where the upper bound is related to energy conservation. The numerical solution of X s = yields the curves in 
Fig l6.5l The variations in a* d {a) with a and d are small. It is interesting to note that not all solutions have finite 
energy: indeed the constraint J dec 2 /(c) < 00 implies for all d and all a the bound a d (ct) > 2 — v. Consequently, for 
inelastic hard spheres [y = 1) the stationary solution with the tail l|6.7|l has a finite energy for all d > 2, whereas in 
the softer Maxwell models (y = 0) the solutions all carry infinite energy. 

Analytic results for the present model are very limited. Near the elastic limit (q small), one finds the exponent 
a* = 2 — x, where x ~ 0(q), by expanding the terms in X s for small x. The result after some calculations reads 
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where — 0.5771... is Euler's constant and ij)(z) — V (z) /T(z) is the logarithmic derivative of the Euler Gamma 
function, the so-called Digamma function. 

VII. CONCLUSIONS AND PERSPECTIVES 

In the present paper we have studied asymptotic properties of scaling or similarity solutions, F(v, t) = 
(i?o(t)) /(u/uo(t)), of the nonlinear Boltzmann equation in spatially uniform systems composed of particles with 
inelastic interactions for large times and large velocities. The large t— and c— scales are relevant because on such 
scales the universal features of the solutions survive, as the velocities, c = v/vo(t), are measured in units of the r.m.s. 
velocity or instantaneous width Vo(t) of the distribution. It follows from DSMC solutions that the v.d.f. F(v,t) 
for large classes of initial distributions, possibly driven by energy sources, evolve into a scaling solution, and after a 
sufficiently long time the combination VQ(t)F(v,i), plotted versus c = v/v (t), can be collapsed for all t on a single 
scaling form f(c) with an overpopulated high e nerg y tail. This observation has been rigorously proven for Maxwell 
models in Ref. [lTj , and for hard spheres in Ref . [23 ■ 

In this article we have focused on the properties of scaling solutions, and in particular on its high energy tail. In 
Section |H] we have introduced the Boltzmann equation for classes of inelastic interactions, corresponding to pseudo- 
repulsive power law potentials, V(r) ~ l/r™, with collision rates scaling like K ~ g v where g is the relative speed of 
the interacting particles. These models embed hard scatterers like inelastic hard spheres (y = 1) and soft scatterers 
like pseudo-Maxwell molecules (y = 0), and even softer ones with v < 0. The energy loss in an inelastic interaction 
is proportional to the inelasticity, (1 — a 2 ), where a is the coefficient of restitution. We also show in Section ITU how 
driving forces, such as white noise or nonlinear friction can be included in the Boltzmann equation. 

Section ITTTI and Appendices lAl and FBI provide the crux of the new method, that enables us to determine the singular 
large— c part h(c) of the scaling form f(c) = 5(c) + h(c). The basic observation is that /(c) = 8(c) is invariant under 
collisions, i.e. I(c\8) = 0. The subsequent linearization in Eqs. i|2.10|l - i|2.11|l of the nonlinear Boltzmann equation 
around this singular stationary solution provides the linearized collision operator A that determines h(c). 

In the Fourier transformation method -which can only be applied to Maxwell models- the corresponding Fourier 
transforms are f(k) = 1 + h(k) and I(k\f), where k is the variable conjugate to c. Then, one has I(k\l) = 0, and it 
leads to the linearized operator I(k\l + h) ~ — A'h(k), where h(k) has a leading small k singularity of 0(k s ), where s 
differs from an even positive integer. This method leads directly to the eigenvalue equations, solved in Refs. 0, 
to determine the power law exponent a in /(c) ~ l/c a+d for Maxwell models. The arguments above, in the reverse 
direction, have generated the essential suggestion to introduce A for general soft sphere models with v 7^ 0, and to 
study its properties. 

Section II VI provides the important link between on the one hand the stability (b > 0) and the marginal stability 
(b = 0) of fixed point solutions of the energy balance equation, and on the other hand the approach of solutions of the 
Boltzmann equation to a scaling form /(c) ~ exp(— /3c b ) with the stretched exponent b > 0, as well as the existence 
of power law solutions /(c) ~ c —a—d—u m t ne i/-model at its stability threshold (b = 0). 

The stretching exponent for the z/-models with white noise driving, discussed in this paper, is b = 1 + vjl. The 
6-values for free cooling and general nonlinear friction were mentioned in the introduction. The existence of power 
law solutions in the corresponding threshold models (6 = 0), as well as the calculation of the power law exponents 
will be published elsewhere ^1 • The application of the general theory to white noise driving for stable soft sphere 
models was presented in section [V] where also the importance of sub-leading corrections was demonstrated in Fia l5.4l 

We have also analyzed an inelastic d-dimensional BGK model [261 l52l | for free cooling and white noise driving. In 
the former case the scaling solution can be solved exactly. The model shows an algebraic tail /(c) ~ c~ a ~ d with an 
exponent a ~ 1/(1 — ce 2 ), where 1 — a 2 is again the fractional loss of energy in a collision. For white noise driving, 
one finds asymptotically /(c) ~ exp(— (3c). So the BGK model exhibits the generic behavior of a Maxwell model. 

In section IVTI the power law exponent a = ad(ct) in the threshold model (b = 1 + vjl = 0) for white noise driving 
is determined from a transcendental equation that can be solved exactly for d = 1 and leads to the power law tail 
/(c) ~ c —a—d-u ^ c -2 ^ ag in us t r ated in Fig l6.ll For general dimensions, the transcendental equation has been studied 
analytically, graphically and numerically and the power law exponent Od(a) for d = 2 and a = 1 is compared with 
the asymptotic result for large d in Fig l6.3l In our new asymptotic method the power law tail for a system driven 
by the ultra- high energy source is simply determined by the eigenvalue equation, Ah s (c) = X s h s (c) = 0, where the 
eigenvalue X s — at s = a*, and where the tail is given by the R-eigenfunction, h a *(c) ~ cr a "~ d ~ v . 

The universal feature of thermodynamic equilibrium in systems with conservative interactions is the Gibbs' state 
in which the v.d.f. is always Maxwellian, exp[— c 2 ]. The conclusion from our analysis tends to be that the only 
generic feature of the present scaling solutions for systems with dissipative interactions is overpopulation of the high 
energy tails relative to Maxwellians. These over-populated tails come in two shapes, either stretched exponentials, 
/(c) ~ exp[— [3c b ] (see Fig l5.1l and Fig l5.3|l . or power law tails, /(c) ~ c~ s with s = a + d + u (see Fig lG.ll and Fig l6.4[l . 
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The stretching exponent b = 1 + \v, and the power law exponent s = a + d + v with a = a^(a) depend linearly on 
the interaction exponent v, and they depend sensitively on the inelasticity and on the driving device used |42|. 

Furthermore, to test the theoretical predictions about tails in the velocity distributions by laboratory experiments, 
there is the additional problem that the fundamental inter-particle interactions, as well as the interactions between 
granular particles and macroscopic driving devices are not really known. 



APPENDIX A: BOLTZMANN COLLISION OPERATOR 



1. Nonlinear operator / 

A convenient representation of the nonlinear collision operator is obtained by changing integration variables w = 
(w ± , ?«||) in l|2.4[) into (uj_,U||) = (w",wjj*), where w** = w*(l/a) according to (|2.3f) . This implies 

Vj* = Vx, V** = -(-£?W|| +PU>||) 
U ± =w"=W ± , U|| = W** = — (pU|| - QTO||), (A.l) 

where p = 1 — q = i(l + a). Consequently <iw = (a/q)du and vn — w\\ — (a/q)(u\\ — v\\). Inserting these results in 
(12. 4|) gives the following representation of the nonlinear collision term, 



I(v\F) = du 



F(u) (A.2) 



with g = u — v. By substituting F(v) = <5(v) and using the explicit form of the collision rate, K = g u \g\\ \ a — g v ~ <T \g\\ |°\ 
it is easy to verify that I{v\S) = 0, at least for v > 0. As we have seen in Section II, negative ^—values correspond 
in the case of elastic interactions to a potential V(r) ~ l/r™ with exponent 0<n<2(d— 1), representing weak 
interactions. To fix the divergence difficulties in K at g = we introduce for v < a small— g cut off in the collision 
frequency, i.e. K oc (g + vo(t)A) v , where A is a positive constant of order unity. All applications in the present 
paper of the cut off K for v < concern NESS's for driven cases, where ^o(oo) is finite and non- vanishing. So vqA 
may be simply replaced by A. The reason for adding the factor vq is solely for theoretical convenience, and ensures 
that the Boltzmann equation after rescaling leads to time-independent scaling equations for models with v < 0. For 
asymptotically large velocities, v = vqC with c> 1, the w— integration in I(v\F) is typically restricted tow < vq, and 
the collision frequency becomes K ~ (v + Uo^) _ ' !/ '|wi|| <T . Consequently, the collision term with the cut off collision 
frequency also satisfies the relation, I(v\8) = 0, because the product K (v±, vu)6(v) is well-behaved, the gain and 
loss term cancel, and the linearized collision operator A in (I2.11|) is well behaved for v < as well. Moreover, for 
asymptotically large— v the collision operator / and A are independent of the cut off. We also note that for v < the 
collision frequency K at large v is vanishing like 1/v'"', and the mechanisms for randomizing the large—, particles 
and redistributing their energy over all particles is essentially lacking. Consequently the tail distribution is expected 
to be heavily overpopulated. 

In our DSMC method for numerically solving the nonlinear Boltzmann equation for negative v— values, a similar 
small velocity cut-off in the collision rate is required as well. 



(3 S = / |an| s = ^ ! = v i ' v f , (A.3) 



2. Evaluation of eigenvalue A s (cr) 

The evaluation of the eigenvalue starts from expression (|3.7fl containing four terms. The first, second and fourth 
one are determined by the angular average of | cos 8\ s , i.e. 

J^ /2 d0(mn0)*- 2 |cos0|'' _ r(^±±)r(f) 
/; /2 d0( S in0)*-2 "r(^)F(i) ! 

where a is a fixed direction. This integral converges for a > — 1. The third term, denoted by \i a \ is more complicated, 
i.e. 

Ai a) = / |c||r|£x+?C||n| s = [ IcHHci+^fr 72 

•/ ri in 

J " /2 d9 (sin Of- 1 [l-z cos 2 9] 3/2 (cos 6) a 
fj 2 d9(sm6) d -i 
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In the last equality we have used the relation c\ + c| = 1 and defined z = 1 — q 2 . Next, we introduced the change of 
variables fi = cf, = cos 2 9 and d[i = —2 cos 9 sin9 d9, to obtain 



Pa 



/ s a + 1 a + d 

2*1 



(d-3)/2/ 



ss/2 



2' 2 



(A.5) 



where B(x, y) — T(x)T{y) /T(x + y) . The integral can be identified as a representation of the hyper-geometric function 



2 F 1 (a,b;c\z) = ^ 



{a)n{b)nz n 
n!(c)„ 



= (l/B(b,c-b)) dxx 13 - 1 (l-x) c - b - 1 (1- zx)- a , 
Jo 

where (a)„ = T(a + n)/T{a) and (j a is given by (|A.3jl . Combining results yields the eigenvalue 

A» = ^{1 + ^,0-2^1 H,^;^!!-? 2 )} -p s /W 



(A.6) 



(A.7) 



3. Linear operator A and its eigenfunctions 

To construct A, defined in (|2.11|l . we start from the representation I(c\8) in (|A.2(1 with F(c) = 5(c) + h(c), and use 
I(c\5) = 0. This yields after some transformations, 



A = A<°> + + A<'\ 



(A. 



where 



A®h(c) 
A^h(c) 
A (b) h(c) 



K(cj_, c\\)h(c) + <5(c) / / duK(ux,u\\)h(u) 



1 



n q 
1 



31 
q 



H 








) 




q 





/ Cll — Mil , 

AT u ± , J ^ ) ft 

P V P 



ui+n 



C|| - qu\\ 



P 



S(c± + nttii 



y y du A(uj_, ii||) h(u)5(c — np r i 



(A.9) 



Inspection of the expressions for A( a ) and A( b ) shows that Ac r = fi r c r+ " , where the eigenvalue \x r is different form A r 
in l|A.7(l . Rather than first calculating [i r explicitly, and then determining r such that [i r — X s , we simply state that 
r = — s — d + 2, i.e. 



Ac- 



Ac c 



-s — d 



(s > 0,c> 0), 



(A.10) 



and verify a posteriori that A s is equal to the expression (|A.7|I . 

Consider first A'°) in l|A.9jl . and substitute h(c) = c~ s ~ d ~ v . By performing similar substitutions as in the steps 
(|A.3jl to \A.h\ we obtain 



aw 



-cr-1 



-Paq" 7 - 1 2*1 



-{s+d+<j)/2 



s - 


\-d-\ 


-a cr + 1 a + d 


z 




2 


' 2 ' 2 


z- 1 



= /8a 2*1 



S (7+1 (7 + d 

2' 2 ' 2 



(A.ll) 
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where z = 1 — q 2 . The last line has been obtained with the help of Eq. (9.131.1) of Ref. [53, and is identical to the 
term containing the function 2 F\ in i|A.7|) . To evaluate A^h(c) we use the relation (|A.9|I . 

*(c-npuy) = _J_*(c-n)^|tt|,|-£), (A.12) 

and integrate out both delta functions with the result, 

A ^Mc) = [ du ± [u"-°h(u)} Hl=c/p . (A.13) 

Setting h(c) = l/c s+d+lf and changing variables = xc/p, we obtain 

^--tS/ptf- (A - 14) 

Here we used the relations (JA.3I1 and a = (s + d + a)/2 in agreement with the last term in (|A.7|I . One similarly verifies 
that A (i) in (|A.8|I for c > yields 

gS+d+v J pS+d+v J c s+a ^s+d v ' 

In summary, we have obtained for c > 0, s > 0, 

Kc- s - d - v = \ s {a)c~ s ~ d and A 1 ^ 8 = A s (cr) c s+!y (A.16) 
and for s = one can verify that, 

A<5(c) = and A f l = (A (ct) = 0). (A.17) 

APPENDIX B: ESTIMATES OF I FOR EXPONENTIALLY BOUNDED FUNCTIONS 

The goal of this appendix is to improve the large— c estimate of the nonlinear collision operator, I(c\f) ~ —f3 a c l/ f, as 
given in the literature [lall8| . The basic idea of the method is to split the v.d.f. /(c) = fo(c) + h(c) into a regular bulk 
part /o(c), say of Gaussian shape, and a small singular part h{c). The bulk part carries the mass of the distribution, 
i.e. J dcfo(c) = 1. Then fo(c) vanishes effectively beyond the thermal range, say c > 3, where the small singular part 
h(c) is living (see DSMC results in Fig l5.2l for d = 1, and Fig l5.3l for d = 2). For c>l the bulk part can be viewed 
in zeroth approximation as a normalized Dirac delta function (5(c), and we represent the v.d.f. as /(c) = 6(c) + h(c). 
Moreover, this caricature of the bulk part, 5(c), makes the nonlinear collision term vanish according to 1)2. 9|l . and 
we can linearize the large— c form of I(c\8 + h) = —Ah(c). Inspection of the contributions A",A^, and A^ in 
(|A.8(l - ljA.9|) in Appendix A shows that the term in A^h, proportional to 5(c), and A^h cx J 5(cj_), are short range, 
and consequently the collision term for c> 1 reduces to its asymptotic form, 

4o(c|5 + /i) = -A 00 /i(c) 
= -K(c ± ,c ]] )h(c) + (l/q)J n K(c x ,c ]] /q)h(\c x +nc ]] /q\) , (B.l) 

where J K(c±, cm) = c v (3 ai as implied by (|A. 15|) . We need an asymptotic estimate of the second term on the r.h.s. 
of 1B.1|) . A( a ) , when acting on exponentially bounded functions h(c) of the form l|3.16[) . 

A^f ~ A^c* exp[-/3c & + (3'c b '} = J(c)c v /(c) (B.2) 

with b > 0. Here the last term denotes the outcome of the calculations. To calculate J(c) we introduce the variables, 
— c^_ = cos 2 8 = and £ = (1 — q 2 )/q 2 . Then the integrand in (|B.1|) is a function if(|cii|) = Kf of |cm| with, 



c 



2 



\c ± + nc {] / q \ = [cl + c 2 /q 2 Y' 2 = c[l + C^/ 2 

K{c ± , H /q) = \c {l /q\ a {cl + c 2 / q 2 }^-°V 2 = q -°c»^ 2 {l + Cri {,J - ay 

/(|cx+nc||/g|) = c^(l + CM) x/2 exp[-/3 C fc (l + C/i) b/2 ] 

, _ /; /2 ^(sin0) rf - 2 g(cosg) _ fidrtl-riV-WyL-WHjy/jl) , 
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where B(x, y) = T(x)T(y)/T(x + y). With the help of l|Rl|l - l|R3|) we obtain, 



^( c ) ^ g+ln l-i i, /' Ml ~ M) (d ~ 3)/ V ( ^ 1)/2 (l + C^- a+x)/2 E^)}, (B.4) 
f +1 B(V>l) Jo 

where E(fi) is defined as, 

E(p) = exp{-/? C b [(l + C/i) b/2 -l]+/5"c b '[(l + CM) b ' /2 -l]} 

~ exp[-Xfi + X'fi] (/x^O) (B.5) 

with X = ^(3b(c b and X' = ^[3b'C,c b . For c » 1 the integrand in (|B.4|) vanishes exponentially fast, unless /i = 
cos 2 6> ~ (grazing collisions). Near /i = factors of the form (1 + A/j,) a can be approximated by unity, and 
itself by the second line of (|B.5|I . Changing integration variables X/i — i, and taking the large— c or large— X limit 
yields for the integral, 

.i 

dnn {a - 1)/2 e- Xl " +x ^ = r((<r + l)/2)/X (a+1 ^ 2 {l + 0(1/X)}. (B.6) 



Consequently the large— c behavior of J{c) is, 

'2±±\ / o \ ( cr + 1 )/ 2 



Inserting (3 a from (| A. 3|> and £ = (1 — q 2 )/q 2 we obtain for the coefficient, 



By combining l|B.l|) . (|B.2ll and (|B.7|I . and observing that /(c) ~ /i(c) for c ^> 1, the nonlinear collision term takes the 
asymptotic form, 

IooW) = -ArC^l - /C CT c- a ]/(c), (B.9) 

with a = + 1)- The derivation also shows that for d = 1 the algebraic correction term in l|B.8|) is absent because 
the angular integral J is missing. So for |c| > 1 we have in one dimension, 

A/i(c) = \c\ v {h(c) - g~ w_1 ft(c/?) - ^"-^(c/p)}. (B.10) 

Here g = 1 — p = |(1 — a) and < g < ip < 1. So, as long as a < 1, the contributions from h(c/q) and h(c/p), 
originating from the gain term, are exponentially separated from the loss term in case h(c) ~ exp[— f3c b ] with b > 0. 



ACKNOWLEDGMENTS 



It is a pleasure to dedicate this article to Carlo Cercignani, on the occasion of his 65-th birthday, in honor of his 
great contributions to kinetic theory in general, and to the study of inelastic Maxwell molecules in particular. We 
also acknowledge useful discussions with Ricardo Brito, Eli Ben-Nairn and Paul Krapivsky. M.H.E. is supported by 
Secretaria de Estado de Education y Universidades (Spain) and the research project FIS2004-271. E.T. thanks the 
European Community's Human Potential Program under contract HPRN-CT-2002-00307 (DYGLAGEMEN). 



[1] C.S. Campbell, Annu. Rev. Fluid Mech. 22, 57 (1990). 

[2] N. Sela and I. Goldhirsch, Phys. Fluids 7, 507 (1995); I. Goldhirsch, M.L. Tan and G. Zanetti, J. Scientific Computing 8, 
1 (1993). 

[3] A. Goldshtein and M. Shapiro, J. Fluid Mech. 282, 75 (1995). 

[4] J.J. Brey, J.W. Dufty and A. Santos, J. Stat. Phys. 87, 1051 (1997). 



19 



I. Pagonabarraga, E. Trizac, T.P.C. van Noije and M.H. Ernst, Phys. Rev. E 65, 011303 (2002) and cond-mat/0107570 
Theory of Granular Gas Dynamics, Th. Poschel and N.V. Brilliantov (Eds), Lecture Notes in Physics (Springer- Verlag, 
Berlin, 2003). 

Granular Gases, Th. Poschel and S. Luding (Eds), Lecture Notes in Physics (Springer- Verlag, Berlin, 2001). 

A.V. Bobylev, J. A. Carrillo, and I.M. Gamba, J. Stat. Phys. 98, 743 (2000); erratum: J. Stat. Phys. 103, 1137 (2001). 

J. A. Carrillo, C. Cercignani, and I.M. Gamba, Phys. Rev. E 62, 7700 (2000). C. Cercignani, J. Stat. Phys. 102, 1407 

(2001). 

E. Ben-Nairn and P. Krapivsky, Phys. Rev. E 61, R5 (2000). 

A. Baldassarri, U. Marini Bettolo Marconi, and A. Puglisi, Europhys. Lett. 58, 14 (2002); A. Baldassarri, U. Marini Bettolo 
Marconi, and A. Puglisi, Math. Mod. Meth. Appl. S. 12, 965 (2002). 
A. Baldassarri, U. Marini Bettolo Marconi, and A. P uglisi, See |(|. 
E. Ben-Nairn and P. Krapivsky, See |(j and cond-mat/0301238 

P.L. Krapivsky and E. Ben-Nairn, J. Phys. A: Math Gen. 35, L147 (2002); E. Ben-Nairn and P. Krapivsky, Phys. Rev. E 
66, 11309 (2002). 

M.H. Ernst and R. Brito, Europhys. Lett. 58, 182 (2002); M.H. Ernst and R. Brito, J. Stat. Phys. 109, 407 (2002) and 
|cond-mat/0111093| 

M.H. Ernst and R. Brito, Phys. Rev. E 65. 040301(R) (2002). 

A.V. Bobylev and C. Cercignani, J. Stat. Phys. 110, 333 (2003); A.V. Bobylev, C. Cercignani and G. Toscani, J. Stat. 
Phys. Ill, 403 (2003). 

T.P.C. van Noije and M.H. Ernst, Granular Matter 1, 57 (1998), and|cond-m at/9803042 
J.M. Montanero and A. Santos, Granular Matter 2, 53 (2000) and cond-mat/0002323 
S.E. Esipov and T. Poschel, J. Stat. Phys. 86, 1385 (1997). 

A. Barrat, T. Biben, Z. Racz, E. Trizac, and F. van Wijland, J. Phys. A: Math. Gen. 35, 463 (2002) and cond-mat/01 10345 
J.J. Brey, M.J. Ruiz-Montero and D. Cubero, Phys. Rev. E 54 3664 (1996); J.J. Brey, D. Cubero and M.J. Ruiz-Montero, 

Phys. Rev. E 59, 1256 (1999). 

I.M. Gamba, V. Panferov and C. Villani, Commun. Math. Phys. 246, 503-541 (2004), and math-ph/0306014 

A.V. Bobylev, I.M. Gamba and V. Panferov, J. Stat. Phys.116, 1651 (2004). 

I.M. Gamba, S. Rjasanow and W. Wagner, to appear in Math. Comp. Model, 2005. 

M.H. Ernst and R. Brito, See @ and cond-mat/0304608 

T.P.C. van Noije, M.H. Ernst, E. Trizac and I. Pagonabarraga, Phys. Rev. E 59, 4326 (1999). 

D.R.M. Williams and F.C. MacKintosh, Phys. Rev. E 54, R9 (1996). 

R. Soto, M. Mareschal and M. Malek Mansour, Phys. Rev. E 62, 3836 (2000). 

N.I. Chernov, G.I. Eyink, J.L. Lebowitz, and Ya. G. Sinai, Comm. Math. Physics, 154, 569 (1993). 

D. J. Evans and CP. Morriss, Statistical Mechanics of Non-equilibrium Liquids, (Academic Press, London 1990). 

E. Ben-Nairn and J. Machta, Phys. Rev. Lett. 94, 138001 (2005). 

Th. Biben, Ph. A. Martin and J. Piasecki, Physica A 310, 308 (2002). 

E. Barkay, J. Stat. Phys. 115, 1537 (2004). 

F. Rouyer and N. Menon, Phys. Rev. Lett. 85, 3676 (2000); D.L. Blair and A. Kudrolli, Phys. Rev. E 64, 050301(R) 
(2001); C. Bizon, M.D. Shattuck, J.B. Swift, and H.L. Swinney, Phys. Rev. E 6, 4340 (1999); J.S. Olafsen and J.S. Urbach, 
Phys. Rev. Lett. 81, 4369 (1998); Phys. Rev. E 60, R2468 (1999); W. Losert, D.G.W. Cooper, J. Delour, A. Kudrolli and 
J. P. Gollub, Chaos 9, 682 (1999); I.S. Aranson and J.S. Olafsen, Phys. Rev. E 66, 061302 (2002). 

A. Barrat and E. Trizac, Eur. Phys. J E 11, 99 (2003). 

S. Chapman and T.G. Cowling, The Mathematical Theory of Non-uniform Gases (Cambridge University Press, third 
edition, 1970). 

P. Resibois and M. de Leener, Classical Kinetic Theory of Fluids (Wiley, New York, 1977). 

C. Cercignani, The Boltzmann Equation and its Applications (Springer Verlag, New York, 1988). 

G. E. Uhlenbeck and GW. Ford, Lectures in Statistical Mechanics (Americam Math. Society, Providence, Rhone Island, 
1963). 

Proceedings "Southern Workshop on Granular Materials", December 10-13 (2003), Pucon, Chile: 
www.dfi.uchile.cl/swgm03 (see 'Contributions', M.Ernst [zip]). 
M. Ernst, E. Trizac and A. Barrat, in preparation. 

G. Bird, "Molecular Gas Dynamics" (Oxford University Press, New York, 1976) and "Molecular Gas Dynamics and the 
Direct Simulation of Gas flows" (Clarendon Press, Oxford, 1994). 

H. Waldmann, in: Handbuch der Physik, Vol 12, S. Fliigge (Ed.) (Springer- Verlag, 1958). 
S. Ulam, Adv. Apll. Math. 1, 7(1980). 

A. Santos and M. H. Ernst, Phys. Rev. E 68, 011305-1-11 (2003) cond-mat/0302285 . 
M.H. Ernst, Phys. Rep. 78, 1 (1981). 

A.V. Bobylev and C. Cercignani, J. Stat. Phys. 106, 547 (2002). 

For certain types of driving the regularity of the solutions /(c) has been shown in Ref.|48| and Ref.|23| for inelastic Maxwell 

models and hard spheres respectively. 

P.K. Haff, J. Fluid. Mech. 134, 40 (1983). 

S. McNamara and W.R. Young, Phys. Fluids A 4, 496 (1992). 

J.J. Brey, F. Moreno and J.W. Dufty, Phys. Rev. A 54, 445 (1996). 

I. S. Gradshteyn and I.M. Ryzhik, Table of Inteqrals, Series, and Products, Alan Jeffrey and Daniel Zwillinger (eds.) Sixth 



20 

edition (July 2000). 



21 




FIG. 3.1: Eigenvalue spectrum \ s (a) for s > of the collision operator in the inelastic soft sphere models (a, v), and shown 
for various values of the coefficient of restitution a. The ordinate shows \ 3 (a)/(3cr for d = 2, a = v = 1, which approaches 
1 for s — > oo, and —1 for s — ► + , whereas A (ct) = 0. The bullet is centered at (2,0) showing that a* < 2, where A a * = 
(intersection point with s-axis), and also that s* — > 2 as a — > 1 (energy conservation in elastic case). 




FIG. 5.1: The figure corresponds to a stable non equilibrium steady state at a = and v — { — 1,0,1,2,3} for white 
noise-driven 1-D systems. The DSMC data show that In /(c) is linear in c b with 6=1 + i//2 > 0. The insert with In /(c) 
versus c shows the highest over-population at^ = — lorb=i (outercurve) , whereas the curve for b = v — 2 is a Gaussian. 
Underpopulation occurs for v > 2. 
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FIG. 5.2: Case of a white-noise driven 1-D system, with a — and a = v — 1. The dotted line displays /(c) vs. f3c b while 
the dashed line corresponds to c" / ' 4 /(c) versus (3c b (6=1 + v/2). Taking into account the sub-leading correction c v ^ 4 yields a 
better agreement with the theoretical result, however the figure clearly shows that this correction is very small. 




FIG. 5.3: Case of a white-noise driven 2-D soft sphere model a = v = —1/2. The DSMC data for In /(c) vs c b with 
6=1 + ^/2 = 3/4, confirm the stretched Gaussian tails exp [-/3c 6 ], where increases with a. The insert with In /(c) vs c 
shows the overpopulation in the tail. 
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FIG. 5.4: Illustration of the relevance of the sub-leading correction (white noise, d = 2, a — v — —1/2). Graph a) corresponds 
to a = 0, and graph b) to a = 1/2. The velocity distribution function is shown as a function of c b , together with two asymptotic 
expressions. The first one is the first order prediction exp(— [3c b ), with /3 calculated by DSMC and 6=1 + ^/2 = 3/4 and yields 
a very bad agreement. Inclusion of the sub-leading correction exp(— (3c b + (3'c b ) (still plotted vs. c b ) with calculated by 
DSMC, and 6' = 9/16, gives a very good agreement. Very striking is the fact that in such a plot, /(c) vs c b produces a linear 
high energy tail (in spite of the importance of the sub-leading correction), which would then be well fitted with an effective 
value of (3: /(c) ~ exp(— /3 e ff c b ). As shown here, such an effective value is markedly different from the true (3 in 15. 51 . which 
indicates that any fitting procedure, aiming at computing /3, is doomed to fail. 
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FIG. 6.1: WN-driven system at stability threshold (6 = 1 + v/2=Q) in 1-D. The DSMC data show a power law tail in the 1-D 
soft sphere model (y — —2). It is compared with the analytic prediction /(c) ~ l/c a+d+I/ ~ 1/c 2 (a = 3), shown by the dashed 
line. 




S 



FIG. 6.2: Graphic solution of JOJ for the marginally stable WN-driven 2-D soft sphere model (a = 1, v = — 2) at a = 
(upper two curves) and at a = 0.9 (lower two curves). The solid curves represent y s — A s (l)//3i and the dashed curves 
J/s = A2(l)s 2 /(8/3i) = pqs 2 /6. The largest intersection point, indicated by an open circle, gives the exponent a = a,d(a) of the 
high energy tail / ~ l/c a+i+v . 
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FIG. 6.3: Exponent a<j(a)/d vs a of the power law tail in the 2-D WN-driven soft sphere model (a = 1, v = — 2). The solid 
line shows the value for d — 2 obtained from (16.41 (see also Fig l6.2l for the graphical solution), and the dashed line shows the 
large d— result x+ = a/d, obtained from 16.61 . 




FIG. 6.4: 2-D counter part of Fig l6.ll Measured power law tail of the velocity distribution in the 2-D WN-driven soft sphere 
model (cr = 1, v = —2) at various a. The dashed lines correspond to the predicted power laws with exponents a,d(a) given in 
FigO 
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FIG. 6.5: Solution of A a (cr) = a) as a function of the restitution coefficient a and b) as a function of dimension d (lower 
graph). The root of A a (cr) = is denoted by a* = a* d (a). 



